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Abstract 

We present a numerical study of turbulence in Bose-Einstein condensates within the 3D Gross- 
Pitaevskii equation. We concentrate on the direct energy cascade in forced-dissipated systems. We 
show that behavior of the system is very sensitive to the properties of the model at the scales greater 
than the forcing scale, and we identify three different universal regimes: (1) a non-stationary regime 
with condensation and transition from a four-wave to a three-wave interaction process when the 
largest scales are not dissipated, (2) a steady weak wave turbulence regime when largest scales are 
dissipated with a friction-type dissipation, (3) a state with a scale-by-scale balance of the linear 
and the nonlinear timescales when the large-scale dissipation is a hypo-viscosity. 
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Experimental discovery of Bose- Einstein condensates (BEC) in 1995 [H |2], some sixty 
years after their theoretical prediction [31 H], sparked a renewed interest in this subject. 
Besides the obvious importance of such systems for fundamental physics, BEC experiments 
provide an excellent opportunity to build and study new nonlinear dynamical systems fabri- 
cated with high degree of control and flexibility supplied by optical means. For the nonlinear 
science and applied mathematics such an opportunity is extremely valuable, because it allows 
to implement and test dynamical and statistical regimes previously predicted theoretically 
and to gain insights about new ones for which the theory is yet to be developed. This 
is because BEC can be described by one of the most important and universal PDE's, the 
nonlinear Schrodinger equation, called in this case Gross-Pitaevskii equation (CPE) [5]: 

where ip is the order parameter indicating the condensate wave function, and T> represent 
possible external forcings and dissipation mechanisms. In general, when JF = and V = 0, 
GPE conserves total energy and particles 

H = J^\Vij\^d^ + J^\^\^d^ = HLiN + HNL, (2a) 

N = l^l^r^x. (2b) 

As GPE model describes a Bose gas at very low temperature, it has been used to study the 
formation of a condensate in [HI [71 [8] . Moreover GPE can be mapped, using the Madelung 
transformation, to the Euler equation for ideal fluid flows with the extra quantum pressure 
term. This is why many concepts arising from the fluid dynamics have been discussed and 
studied with GPE, for example vortices and their reconnection [9j. It was also suggested 
that this model allows statistical motions similar to classical fluid turbulence and a number 
of papers was devoted to flnding Kolmogorov spectrum in such GPE turbulence [101 CH fT2l 

On the other hand, GPE solutions also include dispersive waves which may be involved 
in nonlinear interactions, and an approach known as weak wave turbulence (WWT) can be 
made for GPE. Generally, WWT describes statistics on large ensembles of weakly nonlinear 
waves in different applications, i.e. water waves or waves in plasmas [I5]. Such waves interact 
with each other in a resonant way, e.g. in triads ar quartets, thereby transferring energy 
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(or/and any other invariants) through the scale space forming turbulent cascades similar to 
the classical Kolmogorov cascade in hydrodynamic turbulence. One remarkable property of 
WWT is that, in contrast to hydrodynamic turbulence, power law spectra corresponding 
to such cascades, known as Kolmogorov-Zakharov (KZ) spectra, have been found as exact 
stationary solutions of the corresponding wave kinetic equation [T5] . 

WWT for GPE turbulence was developed in [6l [16] ; the following wave kinetic equation 
was derived 

9ni , f / 1 1 1 1 \ . , 

An / nin2n^n4, \ x (3j 



dt J \ni n2 J 

where = (L/27r)'^(^(k„ t)V^*(ki, t)) is the wave- action spectrum averaged over many re- 
alisations (here L is the size of the bounding box and d is the dimension of the space), 
uji = kf is the wave frequency, and ki = |kj|. Equation ^ conserves the total wave-action 
N = J riidki and the total energy E = J uiUidki which correspond to the GPE invari- 
ants (2a) and (2b) respectively. Besides thermodynamic solutions, equation ^ has two 



non-equilibrium stationary isotropic solutions of the form of n{k) ~ corresponding to 
constant fluxes of energy or wave-action (KZ spectra). In 3D the direct energy cascade 
spectrum has a = 3 and the inverse wave-action cascade has a = 7/3. We will present our 
results in terms of the ID wave-action spectrum n^^{k) = A7Tk'^n{k), i.e. after integration 
over the solid angle. For such a spectrum the WWT prediction for the direct cascade is 
n^^{k) ~ k^^. Note that in hydrodynamic turbulence the results are usually discussed in 
terms of the ID energy spectrum E^^{k) (e.g. Kolmogorov E^^{k) ~ k~^^^). For WWT we 
have the relation E^^{k) = uj{k)n^^{k), which for GPE means E^^{k) ~ k~°'~^^. Further, it 
was predicted in P that in presence of a condensate (due to the inverse wave-action cascade) 
the four-wave resonant interaction will eventually be replaced by a three-wave process with 
an acoustic-type KZ spectrum. 

Previously, there has been a number of numerical simulations of tubulence in 2D GPE 
case and comparisons with the WWT predictions [S], [H]- For the 3D case, we are aware of a 
number of simulations of GPE in freely decaying case, [13], or for the unforced, undamped 
simulation where the initial condition relaxes to the thermodynamic solution, [7] . As far a we 
know no steady state (in the sense of cascade) has ever been reached in any simulation, and 
no direct comparison with the WWT predictions has ever been attempted. The purpose 



of the present work is to revisit the problem of 3D GPE turbulence in the direct energy 
cascade range using the numerical simulations. Our goal will be to find the spectrum and 
to see if (and when) it agrees with the WWT prediction. In cases when the numerics yield 
a spectrum which is different from WWT, we will aim to establish the physical mechanisms 
behind the formation of this state. 

Our numerical domain is a cube with uniform mesh of 256^ points and periodic boundary 
conditions. We integrate equation ([T| by a standard split step method. In order to observe 
the cascade, energy and wave-action are injected directly in Fourier space at wave numbers 
G [9AA;, lOAA;] by a forcing term JF = — z/qc*'^^'^-' with randomly distributed in /c-space and 
time and /o defining the forcing coefficient. To absorb energy at high wave numbers and 
prevent accumulation or thermalization, a dissipative hyper-viscous term V = ii>h{—'V'^)"i/j 
is included in Q, with Uh = 2 x 10~^ and n = 8. 

If GPE is integrated with a forcing at low wave numbers and a dissipation at only high 
wave numbers, it will never reach a steady solution. This is because it admits also the inverse 
cascade of wave-action which will start feeding wave number k = and its close vicinity, 
building a strong condensate Cq = |'?/'(0,t)| that changes the type of interactions from four 
to three- wave resonances [B] . This become clear by looking at figure [l] where we show the 
spectrum at two stages. In the first one, (a), the spectrum exhibits at wave numbers larger 
than forcing a power low close to k~^, consistently with the WWT prediction. As the 
simulation evolves, the condensate grows and the spectrum starts to deviate from the pure 
— 1 scaling: a set of well-defined peaks appears in the spectrum, curve (b). These peaks 
are probably the result of a three-wave interaction similar to ones previously observed for a 
three- wave system in [T7j. In the inset we present the condensate as a function of time: the 
dots (a) and (b) correspond to the instant of times at which the two spectra are computed. 

Regime where the condensate is prevalent was theoretically considered in [IB]- In 
this case, the wave field in ^ can be decomposed as il){^,t) = c{t) + 0(x, t), where 
represents small fluctuations, i.e. <^ c. The condensate part evolves as c{t) = cqc*'"'*, with 
Po = 1^(0, Linearizing the system the new dispersion relation Lj{k) = po± k^k'^ -|- 2po 
can be found, known in the literature as Bogoliubov dispersion. In figure [2] we present 
the numerical evaluation of the dispersion relation taken at the final stage of simulation, 
case (b); results show excellent agreement with the theoretical Bogoliubov curve. For very 
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FIG. 1: Spectrum n^^{k,t) at t = 1 x 10^ (a) and at t = 6 x 10^ (b). In this simulation the 
dissipation is only hyper-viscosity. Slopes k^^ and /c^^/^ are also plotted. Inset: cq = |'i/'(0, t)| as 
a function of time. 




FIG. 2: Dispersion relation in the simulation with only hyper-viscosity evaluated at final time, see 
figure[T]case (b). The Bogoliubov dispersion curve (only positive branch) is superposed with white 
dashed line. Inset: zoom in the zone of low /c's to observe the condensate (horizontal branch). 

strong condensate, po ^ k"^, the Bogoliubov waves become acoustic and u!{k) = ky/2p^ (in a 
reference frame rotating with the condensate speed Pq). Such acoustic WWT was considered 
in [19] and the respective KZ spectrum, E^^(k) ~ k~^^'^, is very close to Kolmogorov —5/3. 
To make comparison with our results we have to take into account that, in this regime, 
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18] . By looking at the late time spectrum (b) in figure [ij we see that our 
results are consistent with the Zakharov-Sagdeev prediction k'^^"^ for the three-wave acoustic 
turbulence, although the scaling range is tiny and not well developed because of the long 
transient range with peaks [17J. 

Note that the condensate keeps growing in this simulation. In order to avoid such growth, 
a dissipation can be included at wave numbers lower that the ones corresponding to forcing. 
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FIG. 3: Spectrum n^^{k,t) at final stage of simulation in the presence of the friction term. The 
prediction of WWT is also shown. Inset: cq = |V'(0,i)| as a function of time. 



Different options are available. Firstly, we will use a friction-type dissipation which takes the 
form, in Fourier space, V = i^6{k* — k)ilj, where 9 is the Heaviside step function, k* = 9Ak 
corresponds to lowest wavenumber forced and is a friction coefficient which has been set 
to /i = 1 X 10~^. We present our stationary state solution in figure [sj The resulting spectral 
slope is consistent with the prediction of the WWT theory. The growth of the condensate 
is now stopped by friction, as shown in the inset, and transition from the four- wave to a 
three-wave regime is prevented. 

Another common way of damping the low wave numbers consists, in analogy to what is 
done at high wave numbers, in including a hypo-viscosity term T> = iui{'V~'^)"^ip in equation 
Q and suppressing the condensate in Fourier space (mode k = 0). In our simulations, we 
have chosen z/; = 1 x 10~^^ and m = 8. In figure |4] we show the stationary states achieved 
with this new damping term for different forcing coefficient /q. The observed spectrum is 
clearly much steeper than the WWT prediction and it is reasonably fitted by a power law 
for forcing /q in a wide range (two orders of magnitude). It seems that the direct energy 
cascade is strongly influenced by the accumulation of wave-action at wave numbers below 
the forcing. In other words, a sharp dissipative term at low wave numbers can cause an 
infrared bottleneck effect. Similar behavior (steeper spectrum) has been observed recently in 
numerical simulations for water waves [20] . 

To understand these results we try to catch the level of nonlinearity in the system by 
considering the ratio rj = H]^^/ Hli^. Note that integral quantities are not always relevant 
because we are interested at rj in the inertial range and both energies may be strongly 
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FIG. 4: Wave-action n^^{k,tf) spectrum at final stage of simulation with hypo-viscosity for dif- 
ferent forcing coefficient: /o = 0.05 (A), /o = 0.1 (B), fo = 0.5 (C), /o = 1.0 (D), /o = 2 (E), 
/o = 3 (F) . A fc^^ slope is also plotted. 
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FIG. 5: Energy ratio t] = Hj^l/ Huj^ evaluated at steady non-equilibrium state with the hypo- 
viscous dissipation for different forcing coefficient /o (see figure |4] and its label) . 

influenced by what happens, for example, in the forcing or in the low wave number region. 
In the case where WWT prediction are confirmed (figure [s]), r] ^ 1.06; apparently WWT 
condition is not valid in this case but probably most of the nonlinear energy, in Fourier 
space, is stacked at low wave-numbers and so, in the inertial range, the nonlinearity remains 
weak. It is instructive to look now at rj in simulations with hypo-viscosity that give the 
slope. As we can see in figure [5} even by increasing the forcing coefficient fo by two order 
of magnitude, the ratio rj remains of order one. In those cases it is reasonable to think that 
the infrared bottleneck accumulation lead to the growth of the nonlinear terms until they 
become comparable to the linear ones in the inertial range. 
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This observations lead to a "critical balance" (CB) conjecture that the systems saturates 
in a state where the linear and the nonlinear timescales are balanced on a scale-by-scale 
basis. The name CB is borrowed from MHD turbulence were it was originally proposed 
in [21]. Even though not called by this name, CB-like ideas were put forward in the past 
for several other physical wave systems, notably the water surface gravity waves where the 
CB condition leads to the famous Phillips spectrum. Indeed, in the Phillips spectrum the 
wave steepness is saturated by the wave breaking process when a fluid particle cannot stay 
attached to the water surface because its downward acceleration becomes equal to the gravity 
constant, which occurs when the nonlinear time scale becomes of order of the linear one. 

Now we propose a similar CB state may form in CPE turbulence with hypo-viscous 
dissipation (or another kind of low-k dissipation which is sharp enough to lead to the infrared 
bottleneck). Namely, when the low-k range is over-dissipated by strong hypo- viscosity, the 
inverse cascade tendency tends to accumulate the spectrum at low k's until the critical 
balance is reached and the spectrum is saturated. When the size of the nonlinear term, 
locally in Fourier space, becomes of the same order as the linear, which is the critical balance 
condition, the inverse cascade is arrested and the further (infrared) bottleneck accumulation 
is halted. One could also qualitatively view this as a set of nonlinear coherent structures, in 
this case solitons or/and vortices, whose amplitude is limited by the linear dispersion (i.e. 
stronger solitons would break into the weaker ones and incoherent waves). 

We now present an estimate the CB spectrum in the GPE model. Equating the linear 
and the nonlinear terms for equation ([!]) written in Fourier space, we have 

fc'l^l ~ (4) 

which gives for the ID wave action spectrum 

ni^(A;) = Ank^nik) = 47rA;2(L/27r)3|4f ~ k'^. (5) 

Note that in evaluating the nonlinear term in ^ we replaced each integration over k with 
k^ which implies that the modes with wave numbers k are correlated with other fc-modes 
with wavenumbers of the similar values ~ k. This is consistent with our assumption that 
the nonlinearity is scale-by-scale of the same order as the linear term. The k~^ prediction 
in equation ([s]) is consistent with our numerical simulations shown in figure |4| 

Concluding, we have performed numerical simulations of the 3D GPE with forcing and 
dissipation. The direct energy cascade range is strongly influenced by the second conserved 
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quantity, the wave-action N, which has an inverse cascade tendency; results depend on 
how the low /c's are damped. We have observed three different types of universal behavior 
roughly corresponding to situations where the largest scales are either non-dissipative, or 
damped by an efficient (e.g. friction-type) dissipation, or or damped by an inefficient (e.g. 
hypo-viscosity) dissipation. In the first case turbulence is not steady: initial direct energy 
cascade, with a spectrum in good agreement with predictions of the WWT theory, is followed 
by condensation at the largest scales and a transition from a four to a three- wave interactions 
with a clearly Bogoliubov dispersion relation characteristic to this regime. In the second 
case, the wave-action cascade is effectively absorbed so that there is no condensation, and we 
observe a steady state spectrum which is in good agreement with the WWT theory. In the 
third case, the dissipation is not so efficient and an infrared bottleneck forms in the spectrum. 
In this regime we observe a robust steady state spectrum which could be explained by a 
phenomenological "critical balance" proposition where the linear and nonlinear timescales 
are balanced on the scale- by-scale basis. 
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